% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function [f,M] = time_spectra(t,y)
% [f,M] = time_spectra(t,y)
%
% Spectral analysis of the time series (t,y). It is assumed that t is
% equally spaced, so that dt=t(2)-t(1).
%
% f: array of the discrete frequencies [1/s]
% M: amplitude of the spectrum

 dtv = t(2:end)-t(1:end-1);
 if(min(dtv)~=max(dtv))
   warning(['The time step is not uniform: max dt ', ...
     num2str(max(dtv)),', min dt ',num2str(min(dtv))]);
   tint = linspace(t(1),t(end),4*length(t));
   y = interp1(t,y,tint,'linear');
   t = tint;
 end

 dt = t(2)-t(1);
 N = length(t);
 DT = N*dt;
 % Make sure k has the same shape of t
 k = t;
 for j=1:N
   k(j) = j-1;
 end

 % Output frequencies (including the effect of fftshift)
 f = 1/DT * ( k - floor(N/2) );

 % Fourier transform of the time series
 hat_y = DT/(sqrt(2*pi)*N) * fftshift(fft(y));

 % Amplitude of the Fourier coefficients
 M = abs(hat_y);

return
